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Abstract: 

Many problems in early vision are ill posed 1 . Edge detection is a typical 
example. This paper applies regularization techniques to the problem of edge 
detection. We derive an optimal filter for edge detection with a size controlled 
by the regularization parameter A and compare it to the Gaussian filter. A 
formula relating the signal-to-noise ratio to the parameter A is derived from 
regularization analysis for the case of small values of A. We also discuss 
the method of Generalized Cross Validation for obtaining the optimal filter 
scale. Finally, we use our framework to explain two perceptual phenomena: 
coarsely quantized images becoming recognizable by either blurring or adding 
noise. 
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1 Introduction 

Edge detection can be considered as a problem of numerical differentiation 
which can be shown to be an ill-posed problem 3 . One approach to solve such 
an ill-posed problem is to regularize it. Standard regularization techniques 
suggest the use of Gaussian-like filters before differentiation 2,4 . In this paper, 
we address the issue of how to estimate the optimal scale of the filter, that is, 
the amount of smoothing required by the given image data. More precisely, 
we will: 

1. Obtain a regularization filter for edge detection and show that it can 
be seen as a generalization of the Poggio, Voorhees and Yuille filter 
(1985) (PVY filter) 2 . 

2. Understand the role of the parameter A (a problem mentioned by 
T.Poggio and V.Torre 1984 3 ) and find out how to calculate it. We 
will show that the optimal value of A, which controls the size of the 
filter, is related to the signal-to-noise ratio. 

3. Compare the optimal filter with the Gaussian filter and PVY filter. We 
show how stable cr (the size of the gaussian filter) is with respect to A. 
This also corresponds to the stability of the width of the optimal filter 
with respect to the parameter A. 

4. Show that the perceptual phenomena of improved recognizability of 
coarsely quantized images by either adding noise or blurring can be 
understood in terms of scale changes in the edge detector’s filter. 

5. Apply this filter to edge detection. Given the image data, we use 
Canny’s 7 edge detector but using the scale provided by the regulariza¬ 
tion analysis. More precisely, the noise in the image is first computed, 
then the signal-to-noise ratio is obtained and finally the A parameter 
is used to set the scale of the filter used in computing the edges. 

6. Compare the standard regularization techniques with the Generalized 
Cross Validation method for finding the optimal A. 

7. Propose biological mechanisms for selecting the optimal scale. 
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We now give an outline of the paper. Chapter 2 is a brief summary of 
results relating to ill-posed problems, regularization analysis and the varia¬ 
tional method. In Chapter 3 we derive the optimal filter. We show that this 
filter gives the solution of the ill-posed problem of finding the exact signal 
from noisy data. In Chapter 4 we compare the optimal filter, Gaussian filter 
and PVY filter. In Chapter 5 we analyze the role of the parameter A. In 
Chapter 6 we present an implementation of the edge detection results from 
Chapter 4 for computational vision. In Chapter 7 we give a possible expla¬ 
nation for the phenomena of improved recognizability of coarsely quantized 
image when either noise is added or blurring. In Chapter 8 we compare 
the method of Generalized Cross Validation (GCV) for finding A with the 
standard regularization analysis. 


2 Framework for Edge Detection 

Regularization Techniques for Ill-posed Problems 
A problem 

Az = u (2.1) 

for which the class S of solutions z, given A and u , is not compact (changes 
on the right-hand side of the equation can take u outside the set A5) is 
called ill-posed. The approach suggested by Tikhonov to deal with ill-posed 
problems is to construct approximate solutions of equation (2.1) that are 
stable under small changes in the data u. 

If the right-hand side of equation (2.1) is known only approximately, we 
have u(x,y) = ux(x,y) + v(x,y) , where uj(x,y) is the true solution and 
v(x,y) is noise. Then, 


where 


z(u;,i/) 


u(u>, v) 
k(w,v) 


Zt( W, v) + 




Az - j I K(x - £,y - r)z(£,r)d£ dr (2.2) 

and k(u is the Fourier transform of K(x,y). It would seem natural to 
take the solution of equation (2.1) as being 
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since uj(u;,i/) = &(a;, i/)2 r (a), i/). However, this function may not exist since 
the last integral may diverge. Furthermore, even if this ratio does have an 
inverse Fourier transform, the deviation from zero (in the C- or Z 2 -metric) 
can be arbitrarily large, and thus, we cannot think of the exact solution of 
equation (2.1) as an approximate solution of the equation with approximate 
right-hand side. 

Finding edges in an image is in general an ill-posed problem 4 , since it 
involves taking an appropriate derivative of noisy data (notice that we do 
not specify which derivative operator should be used: it may be a directional 
derivative 4 or any other desirable differential operator). The differentiation 
of the function u(x) is ill-posed, since it can be seen as a solution of equation 
2.1 for the operator A of the form 


/ *(£)<# = / Hx ~ = u ( x ) 

J — OO J — oo 

where h(x) is the step function. As described by Rheinsch 9 and by Poggio and 
Torre , this problem can be regularized by smoothing the data before taking 
derivatives. The idea is to consider the regularized solution z of equation 
(2.1), with A being the imaging operator, such that z is sufficiently well- 
behaved for numerical differentiation and as close as possible to the true 
ddtd* 

Tikhonov 1 proves that, for the case of one dimensional image data, to 
approximate a solution of equation (2.1) one takes the solution of a different 
problem, the one of minimizing the functional given by the following equa¬ 
tion, that is close to the original problem for small values of the error in the 
data: 


M x [z,u}= f [ [Az - u] 2 dxdy + Aft[z] (2.3) 

J — oo J —oo 7 

where u(x , y ) is the image data, A is the regularization parameter, A is in 
our special case a convolution operator and fl[z] is the stabilizing operator. 
We will be considering the special case where 
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foo foo \ / ffi ffi \ P/2 

nw = /./-J(a? + aiO z{x ’ y) 


dx dy. 


The operator fi[z] forces to smooth the data and the motivation to choose 
this particular stabilizing operator is to guarantee a desired smoothness in 
the solution. The order of the derivative in ft[z], controlled by the parameter 
p, should be high enough to ensure the appropriate degree of differentiability 
in z required by the derivative operations that has to be performed next. 

It should be pointed out that, whereas the original problem (2.1) does not 
have the property of stability, the problem of minimizing the functional 
M x [z, u ] is stable under small changes in the right-hand side u. This stability 
has been attained by narrowing the class of possible solutions by introducing 
the stabilizing functional fi[z]. 


3 The Optimal Filter 

* 

The Fourier transform of (g- + |^) p/2 z is M(u>, i/)z(w, i/), where = 

(w 2 + i/ 2 ) p . Using Parseval’s theorem we can rewrite equation (2.3) as 

+ A /_ J (w 2 + i ' 2 ) p z(u, v)z(-u, - v)dudv j. 

The associated Euler-Lagrange equation is 

dM x 

dz(-uj - v) = ^ v ) ~ uv )\ k {-u,-v) + A(w 2 + u 2 ) p z(u}, u) = 0. 

For the special case where p = 2 (stabilize with Laplacian) and the operator 

A is given by (2.2) with k(w,u) = e -|(“ a +^), we obtain the regularized 
solution 


z \(x,y) 



_1_ 

1 + A(a; 2 + i/2) 2g6(a; J +i/ J ) 


z(u>, 


v)e-^*+ v ^dudv 


v)\dw du 


4 




4l -' 80 SO 100 


Optimal Filter 


Figure 1: The optimal filter is plotted in the space domain for different values 
of the regularization and scale parameter A 


that corresponds, in the Fourier domain, to filtering the data with 


f{u,v, A) 


1 

1 + A(w 2 + z/2)2 e 6(^+^)- 


(3.1) 


The term k(u),v) approximates the point spread function of the imaging 
device. For the human eye the A(u/,i/) is described by a Bessel function and 
a Gaussian can be a good approximation. In this case the bandwidth (<r) 
of the point spread function of the pupil is of the order of 20 seconds of arc 
which implies b = 0.5<r 2 = 1.5 10~ 5 per degree square. For CCD cameras 
the bandwidth is less than a pixel, restricting the effect of this term. In our 
experiments we used b — 0.5 per pixel square. 

The case b=0.0 gives the PVY filter. In this case the filter is not smooth 
enough to ensure differentiability with a second order differential operator 
such as the Laplacian and a higher order stabilizer is needed 2 . For 6 > 0 
this problem, however, disappears. The filter is plotted (b = 0.5) in the space 

domain for different values of the parameter A (see figure 1), which controls 
the size of the filter. 


Appendix 2 shows how to obtain the Wiener filter using the variational 
method. We can see that the Wiener filter is slightly different from the 
optimal filter, because of the stabilizer term. 
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Figure 2: The relation between A and <7 - the width of the optimal filter and 

4 -Ktn4)^ 

of the Gaussian, respectively. A = f-re 


4 Filters: Optimal, Gaussian, PVY and Wiener 
4.1 The relation between A , a , a and j 

In order to compare the filters, we have to assume an equivalent criterion for 
the Optimal filter, the Gaussian filter, PVY filter and Wiener filter. A pos¬ 
sible criterion is to choose the same half amplitude in the frequency domain. 

1 ) Gaussian-filter - | = e j then u>„ = ( l ~~ 

2) Optimal filter - A = 4 then Aw* = e~ b “* 

3) PVY-filter - } = then aw* = 1 

4) Wiener-filter - \ = e then 7 = e -6w “ 

Then we plot a graph <7 x A for b=0.5 per pixel square (see figure 2 ) 

We notice here the stability of the parameter <7 (the bandwidth,in pixel 
unity,of a filter) with respect to the parameter A. 

Another criterion for equivalence is to chose filters that have the same 
bandwidth in their second-derivative (see figure 4). This is important when 
computing zero-crossings, since in this case the desired operator of convolu- 
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tion is the second derivative of the filter. The value of the parameter A may 
vary with these different definitions, but the actual bandwidth of the filter 
will not vary significantly. 


5 The Optimal Scale 

In order to find the optimal scale we vary the parameter A to obtain the 
closest solution to the true solution. More precisely (see appendix 3), we 
minimize 


' | e -^+.»> + A(u,J + yp- ,L,dv 

(5.1) 

where E{.} is the expectation value operator, S(u>,i/) is the spectral density 
(power spectrum) of zx(x,y) and is the spectral density of v(a:,y), 

assuming that v(x, y ) is stationary. 

Differentiating equation 5.1 with respect to A we obtain 


\ e b{u 2 +l, 3 ) S(w,v)(u 2 + i / 2 ) 4 ; , _ [°° N{u}v)e- b ^ i+l, 2 \<jj 2 + I / 2 ) 2 

Jo Jo (e-H“ 2 +‘' 2 ) + A(w 2 + Z / 2 ) 2 ) 3 dU; dV ~ J 0 J 0 (- e - 6 (^) + A(u; 2 + ^2)2)3^^ 

(5.2) 

For the discrete case the space is defined in a grid with N x M pix- 
els^and the spatial frequency given by = gt and v, = gj where i = 

-y, ...,0,..., J and j = -y,...,0,..., y. We follow the same calculations 
as in the continuous case and keeping the order of the stabilizer a general 
parameter p (before we set p = 2). Using the symmetry with respect to zero 
of 5(w,i/) andN(u,i/) one obtain similarly to equation 5.2 the equation 


£ £ A e~ b2lr W a+( A) i )s(i,j)(2^) 4 P((l)^ + (1)2)2 p _ 
i=o j=o ( e -b 2 »((i)*+c*)») + A(2)r) I P((i)2 + (^P)p) 3 

£ £ N(i,j)e-^((^-K^ ! )(2^)^((i)^ + (A) 3 )" 

i=o j=o (e- b2 "«A)’+(*)’) + A(2jr) 2 P((i)2 + (± ) 2 )p) 3 
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Within the range of the parameter A there are two cases two consider. 
One for very small values, i.e. when A —» 0, and the other case for A not small. 
The second case is applicable whenever the signal-to-noise ratio is small and 
to find A one has to solve equation 5.2b. The first case is applicable whenever 
the signal-to-noise ratio is very small and the following short cut to find A 
can be used. We assume the asymptotic values 


lim S(u!,w ) 

*oo 


So 

(a ; 2 + u 2 ) c 


c > 0 


lim N(u>, u) 

U),V — M3© ' 


N 0 

(a ; 2 -f- u 2 ) a 


and considering the case with white noise (a = 0) we obtain (see appendix 
4) A as the solution of 


a > 0 



A)*-? = In 


! («(^a)A+2 P +i 


(5.3) 

For the case b = 0.0 (the PVY filter) we get A = (2 p + l/4p — . 

A graph of signal-to-noise ratio versus A for b = 0.5 per pixel square and 
is plotted in figure 3. Thus the optimal size of the filter can be obtained 
directly from an estimate of the signal-to-noise ratio for any given image. 
The equation 5.3, for finding the parameter A, turns out to be sensitive when 
we fix the value of the signal-to-noise ratio and vary the values of b or c. 
However the relevant parameter, namely, the width of the filter is not very 
sensitive to it. In our experiments the width did not vary much from 0.8 
pixels to 1.1 pixels. The parameter c was varied from 2 to 15, the parameter 
p (degree of smoothness) was varied from 2 to 6 and the parameter b (optical 
blurring) was varied from 0 to 2. We have to remember that here we are just 
considering large values of signal-to-noise ratio, typically jf- > 100. 


0 Noise Estimation 

Here we assume the noise to be caused by the deficiency in the optical sys¬ 
tem or by artificially adding a random variable to the image. So in this way 
noise is independent of the recognition task and independent of the scene or 
a particular scale. To estimate the noise we use a technique developed by 
H.Voorhees 10 . First, the gradient of the image is computed. A Rayleigh 
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Figure 3: Signal-to-noise ratio versus X, for b = 0.5. The last curve the 
x-axis is given by the width of the filter u 


distribution is then fit to the histogram of the norm of the gradient and the 
noise parameters are estimated. The signal power is obtained from the vari¬ 
ance of the intensity of the image. With this method the program estimates 
from the image data the corresponding signal-to-noise ratio. This ratio gives 
A from relations such as equation 5.3 or 5.2b. The results are shown in 
figures: 4, 5 and 6. One alternative for noise estimation could be done by 
taking a sequence of images and from the temporal correlation obtain the 
noise parameters. Another one is, assuming the noise to be white, to find 
the average value for large frequencies of the image-spectrum, since typically 
for large frequencies the noise is constant and the signal values close to zero. 

7 Two perceptual phenomena: explanations 
and biological implications 

7.1 Coarse quantized images can be better recog¬ 
nized when noise is added 

We first discuss the perceptual phenomena of improved recognizability of 
coarse quantized images when noise is added 5 . Consider an image with 320 
by 384 pixels and 8 bits. A coarsely quantized version of it is shown in 
figure 5a. The optimal filter for figure 5a, estimated as explained above, 
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Figure 4: Zero-crossings , the Laplacian applied to an image that was 
smoothed using a Gaussian filter and the optimal filter with the same width 
for the second derivative 


turns out to have a small scale of % 1.0 pixel corresponding to very large 
signal-to-noise ratio of SJN 0 = 1000.0. We then apply the Canny’s edge 
detector using a gaussian filter with the same bandwidth predicted by the 
regularization. The edges do not easily reveal a face as we see in figure 5c. 
Gaussian white noise with standard deviation 70 is added to figure 5a (see 
figure 5b), making recognition easier. Estimation of the optimal scale, using 
equation 5.2b since the signal-to-noise ratio is small, gives now a width of 
« 6.0 pixels. The corresponding contours (output from Canny’s detector) 
reveal the face in a much better way. 

These results may shed some light on what the visual system may be doing. 
Harmon and Julesz 6 claim that for the quantized image “high frequencies 
introduced by quantized blocking mask the lower spatial frequencies which 
convey information about the face, preventing recognition”. In our frame¬ 
work two processes determine recognizability of the face. The first process 
consists of the estimation of the signal-to-noise ratio (S 0 /N 0 ). The second 
step is to use S 0 /N 0 to set the optimal A for computing the correspon din g 
“edges”. In the case of the quantized image the ratio S 0 /N 0 is large. A is 
then small, which implies that a large bandwidth channel (in the spatial- 
frequency domain) is selected. The zero crossings for this channels do not 
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Figure 5: a,c) Quantized image and edges with <r - 1.0 pixel. b,d) White 
noise, standard deviation of 10 units, was added to the quantized image. 
Now the edges are computed using a — 6.0 pixels. 
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easily allow face recognition because they mostly capture the box outlines. 
For noisy quantized images the ratio S 0 /N 0 is small and correspondingly A is 
large. This implies a filter with small bandwidth. In this case the small band¬ 
width filter suppresses the noise and, as a side effect, also the high frequency 
outlines of the boxes. 

This explanation is not in contrast with the one given by Canny 7 or by 
the one given by Morrone, Burr and Ross (1983) when they claim “that added 
noise (more high frequencies) destroys the propensity to organize the image 
according to its spurious high-frequency structure,...”, but is more precise. 

7.2 Improved recognizability of coarse quantized im¬ 
ages by blurring 

Blurring coarsely quantized images also improves recognition 8 . The explana¬ 
tion for this second perceptual phenomena is natural in our model. Consider 
the zero crossings of the blurred image. For the Gaussian filter they are 
represented by the solution of 

V 2 (G(<r 2 ;x,y) * (G(<7 i; x,y) * I(x,y))) = 0 

Where I(x,y) is the quantized image, * stands for convolution, G(a;x,y) 
is the Gaussian and the blurring process is ((?(< 7 x; *, y) * I(x,y)). Because 
all of those operations are linear we can first convolve both Gaussians. In 
the Fourier domain this is equivalent to multiplying the Gaussians, and the 
multiplication of two Gaussians is a Gaussian with <r = [(<Tx) 2 + (o^) 2 ] 1 / 2 . So 
computing zero crossing of a blurred image is the same as computing zero 
crossings of the same image but with larger <r (a = [(<n) 2 + ((T 2 ) 2 ] 1/2 ). 
Therefore blurring is equivalent to using an effectively larger filter for edge 
detection. This has the effect of suppressing the spurious high frequency 
edges introduced by coarse quantization. 1 

1 Notice that blurring the quantized noisy image has the effect of increasing the esti¬ 
mated signal-to-noise ratio, thereby reducing <r 2 to a value close to the one obtained for 
the quantized image. 
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8 The Method of Generalized Cross Valida¬ 
tion and Regularization 

When S 0 /N 0 cannot be directly estimated, it is natural to consider the 
method of Generalized Cross Validation (GCV) 8 . The GCV method states 
that the optimal value of A can be obtained by minimizing the functional 
(here in one dimension) 


v(x) = -E 

n r-J 


i=0 


[Az n ,\(tj) - Uj] 2 
(l-a* fe (A)) 2 


^ 2 (A) 


(8.3) 


where Wfe (A) = (1 - a**(A))/(l - ± £] =0 a^A)) and a kk ( A) = JL(Az nA )(t k ). 
From here on we use the Gaussian filter , since would be computationally 
more expensive to use the optimal filter. Equation 8.3 reduces to 


V(a) = 


vhJ 2 i= 


—(i —fc ) 2 

e 2»- 2 z(k ) — z(i) 


The method is computationally expensive but intrinsically parallel. We have 
implemented it on the Connection Machine 2 . We tested this method on dif¬ 
ferent images including the ones in figure 5 with various amounts of noise. We 
notice a consistency of the GCV with the results obtained with the standard 
regularization. The GCV gave allways (see figure 6) a slightly mailer value for 
the parameter A (and <x, the width of the filter). We point out that equation 
5.3 is just applicable when the signal-to-noise ratio is high (higher than 200) 
and the other cases we used equation 5.2b. The results of applying C ann y’s 
edges detector with both estimates of the parameter suggest a slightly better 
performance for the standard regularization method. For large amount of 
noise as in figure 5b we obtained significantly different results. For figure 
5c for example the GCV method gave a = 3.5 pixels as oppose to 6.0 given 
by the standard regularization method. This suggests that since no noise 
estimation is involved in the GCV method whenever the signal-to-noise ratio 
is small the standard regularization is going to give better results. Typically 
slices of the image of size 80 pixels require 20 milliseconds for computing 
V(cr). Using Newton’s method to find the minimum, the algorit hm con¬ 
verges after 10 iterations of V(<7). Therefore the GCV method takes in this 


2 The Connection Machine is a massivelly parallel computer with 65,536 processors 
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case 0.2 seconds to find the optimal a. In figure 6 we show several examples 
where the optimal scale was estimated with the standard regularization and 
GCV methods. 

It is interesting to notice that for the figure 6.a the optimal scale was for 
A = 0.1 (cr = 0.91 pixels), however there are two distinct scales for this image 
that could be considered optimal. The coarser scale was not captured with 
any regularization method. A possible way of obtaining both scales would 
be by reconsidering the definition of noise or signal. For instance if the noise 
was defined to be the medium and high frequency of the spectrum then the 
noise estimation would be higher than the one we used and the coarser scale 
would be optimal. 


Figure 6. Several examples of applying the standard regularization method 
and the GCV method to obtain the optimal scale 



Westminst and Canny’s edges with <7 — 0.86 and (T — 0.66. Cross- 
validation: a = 0.66 pixels, standard-regularization: cr - 0.86 pixels 
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Cloth and Canny’s edges with cr = 0.91 and a = 0.54 
= 0.54 pixels, standard-regularization: cr = 0.91 pixels. 


Cross-validation: 


WALTER 





Walter and Canny’s edges with cr = 1.63 and cr — 1.57. Cross-validation: 
a — 1.57 pixels, standard-regularization: cr — 1.62 pixels. We solved equation 
5.b since So/No ratio smaller than 200. 



Apple and Canny’s edges with cr — 1.2 and <j = 0.66. Cross-validation: 
<7 — 0.66 pixels, standard-regularization: cr = 1.2 pixels. We solved equation 
5.b since So/No ratio smaller than 200. 
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Ericqn and Canny s edges with <j — 6.00 and & = 3.5. Cross-validation: 
® 3.48 pixels, standard-regularization: «r = 6.08 pixels. We solved equation 

5.b since So/No ratio smaller than 200. 


9 Conclusion 

We have derived rigorously the optimal way of filtering images prior to nu¬ 
merical differentiation. We also obtained the precise relation between the 
scale of the filter A and the signal-to-noise ratio of the image. Some bio¬ 
logical implications were also considered. In particular we suggested that 
humans can estimate the signal-to-noise ratio in the image from which the 
scale A is computed. Only channels with the appropriate spatial-frequency 
band are then used, the others being inhibited. In this framework it is possi¬ 
ble to understand the perceptual phenomena of improved recognizability of 
coarsely quantized images when noise is added. When the signal-to-noise ra¬ 
tio is large, the estimated A is small and the associated edges do not provide 
good information for recognition. When S 0 /N 0 is smaller, the estimated A is 
larger: the edges then provide better information for recognizing the face in 
the image. When the signal-to-noise ratio cannot be estimated, it is possible 
to use the method of Cross Validation for estimating the optimal A. 
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D.G.). Thanks to E. Hildreth for the suggestions, to E. Gamble, M. Gen- 
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comments. 
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Appendix 1: Deviation of the regularized solution from the ex¬ 
act one and the stability of the solution 

The regularized solution can be written in the form (see Chapter 3) 


i f 00 r 

^ z ’ v) = (wJ-J- 


1 too ,oo 


u T {u,u)e-* wx+l,y) dudv+ 


(2tt) 2 7-00 7-00 K{uj)K{u) 

r r v (“> ^e-^^du dv 

(2x) 2 7-00 7-00 K(u>)K(u) 


then 


* k(x ’ y) - zr(x ' v) = (5 U{ k£)K(v) ] <*+ 


1 /■<*> /(w, V, A) 


— T 

Jir) 2 7- 


(2ir) 2 7-oo A(a;)A(i/) 
Now one defines 


v(aj, t/)e~^ ux+,/t, ^du> dv 


(l-ol) 


^ - (SvC IK kC)ku ( 1 ... 2 ) 


(2tt) 2 7-00 K{u)K{y) 

1 f 00 f 00 /(cu, v. A) 


1 /-OO re 

iN{t ’ X)= (wJ-*,J- 


(2tt) 2 7-00 7-00 A(w)7£T(i/) 


v(u>, i/)e , ( ux+l/ y)duj dv 


(l.a3) 


so that z\(x,y) — zx(x,y) = As(x,y, A) + A^(x,y, A) One notices that equa¬ 
tion (l.a2) characterizes the influence of the regularization and equation 
(l.a3) characterizes the influence of the noise. 

Under these assumptions, the variance of the random function An is 


= „’(A) = JL/~ 
and A 2 (A) = sup t A 2 s (t, A) 

How stable is the regularized solution; more precisely, how much does the 
regularized solution change with a change in the order of the stabilizer p ? 

It is possible to show that: 

1) The asymptotic ( u>, u —>■ oo ) behavior of As and An are the same for 
stabilizers of the form 


M(uj, u) = (u> 2 +i/ 2 ) p +a p ^ 1 (uj 2 +i/ 2 ) p 1 +...-t-a 0 or M(u, u) = (u> 2 +v 2 ) p 
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2) From formula 5.1 


rri\ \_ 1 f /°° A 2 (w 2 + u 2 ) p S(u;,u) , , y* 0 

T(X,p) - — j / r—-;■■■ w 2 — Tr-7zdudt/ + 

2ir 2 Uo [ L(uj , i/) + A(w 2 + i/ 2 ) p ] 2 7o 

= A 2 (A) + v> 2 (A) 


Z(a>)iV(u;) 


[£(a;, i/) + A(u > 2 + i / 2 ) p ] 2 


da; du 


obviously, A is a function of p. Consequently, T(X,p) = ?J>(p) and the function 
VK.P) has a unique minimum at p = p 0 = c — a, and increases monotonically 
in the interval p > p 0 with a finite limit for p —► oo. Tikhonov 1 has shown 
these results for the one dimensional case. 

Defining the measure of stability for arbitrary p > 0 as being 

rP(p) - ^( Po ) 

V’(Po) 

one has from Tikhonov 1 


0 < V’(p) - Hpo) k 7q 

i>(Po) 1 - 7o 

where 7 0 = 2n + 2 c- 2 a an< ^ n * s order of smoothness of the operator K(u>) 
, i.e., how fast A:(u>) —► 0 when u —> oo . 

For the case where k(uj) = e -6 " 2 then n -y oo. Therefore ti slr *<?<>) _> q, 
which is to say that the function V’(p) is weakly dependent on the order 
p of the regularization. Consequently, one can quite justifiably replace the 
optimal order of regularization p 0 with an arbitrary order p > p 0 . For white 
noise p 0 = c and for p — 2 one will be safe whenever working with images 
that at high frequencies do not fall to zero faster than a; -2 . 

It is clear is that if dealing with images with different asymptotic behavior, 
one can reset p to give a good solution, and the stability for higher p will be 
guaranteed. 

In practice we are interested in the stability of the scale of the filter and 
we notice that the actual width of the filter did not vary more than 0.4 pixels 
when the parameters c, p and b (the point spread function parameter) were 
varied. 


18 



Appendix 2: The connection between the regularization method 
and optimal Wiener filtering. 

The regularized solutions of equation 2.1 is 


*>(*»») 


1 k(-u)k(-t/)u(w, v) i( wx+vv )j 

(2tt) 2 7-00 L(uj, v) + AM(w, v ) 


and can be written in the form of a convolution 


/ OO TOO roo 

/ / L\(t - t)u(t)<It 

-OO J — OO <7—00 

where 


Ta(< - r) = 



fc( uj)k( v) c -H^+vy)j^ 
L(u>, i/) + XM(u>, v) 


Now taking T(XM) (formula 5.1), from the condition that this functional be 
minimized on the set of functions v), one finds by elementary calcula¬ 

tions that the minimum is attained with the function 


M(u>, v ) = M 0 ( a>, v) = 


1 N(u), u) 
X S(u>, u) ‘ 


Therefore the approximate (regularized) solution z op (x,y ) of equation 2.1 
obtained is 


Zo P (x, y ) 



Jc(— u)k(— l/)«(cu, V ) 




which is independent of the parameter A. It coincides with the result of 
applying the optimal Wiener filtering to find z{x , y) from u(x, y) = ut(x, y) + 
v(x,y) . 
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Appendix 3: The expression for E{\z\{t) — z-r{t )| 2 } 

We show the result for a more general class of solutions where the smoothing 
term is given by M(uj, u) = ( u 2 +i / 2 ) p . Notice that L{u), u) = k(—u))k(—v)k{uj)k(v). 
We have 


z\{x,y) - zt(x, y) = 

= J—r r { *(-u>)t(-«)M »,v )+»(<■>■«)] . = 

(2x) 2 £(o,,i/) + AAf(w,i/) n ’ " 

= r r, £(<*>, ■<)+ k(-»)k(-v)v( U ) „ )}e ..-(^) (4jdl/ = 

(2tt) 2 J-ooJ-oo' %i/) + AM(w,i/) n ’ W 

- 1 /°° Z 00 AM(w,y)z r (h>,y) t. , v, , 

(2tt) 2 Aoo J-oo ^ L(u>, v) + \M(u>, v) I(u,v) + AM(w,^ ^ + 

since ut{w, v) = k(u>)k(i/)zT(u>, v). Therefore 

E{\z x (x,y) - z T (x,y)\ 2 } = 

= £ /_j_r r 

\(2w) 2 J-oo J-oo 1 L(oj,t/) + 

(27r) 2 7 -co J — co X(u/, l/') + AM(o/, I/') 1 

_ 1 too too too too \ 2 M(u),v)M(<jJ,v')E{zt(u,v)zt(<jJ,v')} 

(47T 2 ) 2 J-oo J-oo J-oo J-oo [L{uJ, v) + AM(w, £/)][X(u/, V 1 ) + AM(ti/, I/')] 

+ 

[X(u>, i/) + AM(w, i/)][£(u/, i/') + A M(u>', i/')] 
since £J{u(a;, i/)} = _E{w(u/, i/)} = 0. For stationary random processes, 

E{z t (u, v).zt{u>', i/)} = S( u>, i/)8(ut + u>')8(i/ + v') 

E{v(u}, l/)v(u!, j/)} = N(u>, l/)8(u! + (Jj')8{v + v')i 


where 8{u> + u/) is Dirac’s delta function. Performing the integration in 
the last expression with respect to u>' and using the properties of the delta 
function and the fact that M(w ) and L(u) are even functions, one obtains 
the value of the deviation 


1 /*oo /*oo 

E{\z x (x,y)-z T (x,y)\ 2 } = T(XM) = — 


\ 2 M 2 (u>, v)S(u}, v) + L(u>, v) j j 

[%.) + AM( W ,.)] 2 
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Appendix 4: Solving equation 5.2, based on Tikhonov 1 

The argument for using the asymptotic expressions (p is the order of differ¬ 
entiation in the smoothing term) follows immediately from 




n o 

. 


e- b ^ 2+ ^S(u, v)(u 2 + v 2 Y 
(e-W+''*) + \( U} 2 + v 2 )vy du}dv 


since <p 2 ( A) —► oo as A —► 0. On the other hand, for any fixed A > 0, 
the integral y? 2 (A) converges. Therefore, for sufficiently small values of the 
parameter A , the basic contribution to <p is provided by large frequencies a?. 
Then one can replace equation 5.2 by 


foo yoo ^ e - 6 ( u ' 2 +,/ 2 )S 0 (u > 2 + v 2 y~ c , , /•<*> /oo e~ b ^ 2 +,/ 2 )N 0 (u? + v 2 y~ a f f 

Jo Jo (e-o^ 2 ) + + v*yy Jo Jo (e-w+" 2 ) + \(u* + v*yy dudl/ 

(5.26) 

with the change of variables uj = £ cosd and v — £ sinO which imply <f 2 = 
u 2 + v 2 and dudv = £ d£ dO . We can rewrite equation (5.2b) as 


r ^ Luos,,?*- 2 '* 1 ,, r ^M 0£ 2p - 2o+1 
Jo (L aa (o + H p y * Jo (uo + w 

where L as (£) = e~ b ?. The integrals 


(5.36) 


h 


=r 


L.M)e* 


—2a+l 


-A 


-r 

Jo 


L..(oe- 


—2c+l 


iwf)+»pr " Jo (iuo+M ”) 3 

are evaluated by the method of steepest descent and are equal to 




Ii = 


h = 


7T 


er" +1 ia S (6) 

{i..te)+Aer}3 1 vi/; , («i.A)i 

-2a+l 


+ 0(AT)] 


M6) 


7T 


{ia.te)+A«r} 3, v !/;■«., A)i 




where 


/,({,*) = Aft. 


{£.,«) + A?} 3 
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and the double prime denote* the i 
& is a root el the equation 


to (. Here, 


(- 2 i„ + - (*r+u)\e-'l„ + » o . 

Substituting the vnluss found lor h and h into equation (2k3b) end keeping 
only the principal terras, 

A 


A - ^[6<A)]-***- fc 


y 
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